###################################################################################################
## Payson (2018) "Cities in the Statehouse" Replication Code
## 50 state analysis in 2007 and 2012 (Table 1 and Figure 1)
## Created with R version 3.4.1 (2017-06-30) -- "Single Candle"
###################################################################################################

## Set working directory
setwd()

## Load packages 
library(ggplot2); library(plyr); library(plm); library(lmtest); library(effects); 
library(stargazer); library(sandwich); library(RColorBrewer); library(reshape);
library(statar); library(fiftystater); library(maps); library(xtable)
data("fifty_states")


## Functions for tables (Note: tables will not compile without these functions)
pretty <- function(x) {
  paste(prettyNum(x, big.mark = ","))
}

clusters <- function(x) {
  length(fixef(x))
}

N <- function(x) {
  dim(x$model)[1]
}


###################################################################################################

## Read in 50-state data
load("50statelobby.Rdata")


###################################################################################################
## Figure 1: Geographic Distribution of Lobbying Cities (Population Over 1,000)
###################################################################################################

## Subset to 2007 to prevent duplication
year2007 <- subset(dta, year == 2007)

## Get state border coordinates
x = map("state", plot = FALSE)$x
y = map("state", plot = FALSE)$y

plot.dta <- data.frame(x=x, y=y)

## Isolate the longitude, latitude, state name, and lobbying status of cities (ever.lobby = indicates
## if city ever lobbied in either 2007 or 2012)
plot.dta1 <- data.frame(x=year2007$long, y = year2007$lat, 
                        lobby = factor(year2007$ever.lobby), state = tolower(year2007$state))

## Subset to cities that ever lobbied 
plot.dta1 <- subset(plot.dta1, lobby == 1)

## Figure 1: Geographic Distribution of Lobbying Cities
ggplot(plot.dta, aes(x=x, y=y)) +
  geom_path(alpha = .5) +
  geom_point(data = plot.dta1, aes(x=x, y=y), 
             alpha = .5, size = 1) +
  theme_minimal(base_family = "serif") +  
  theme(plot.background = element_blank(),
        legend.position = "top",
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())



###################################################################################################
## Table 1: City Lobbying and State Revenue Transfers Across the U.S., 2007 and 2012
###################################################################################################

load("50statelobby.Rdata")

## Format 50-state panel data using plm package
## cog.id = individual Census of Government city ID
pdta <- pdata.frame(dta, index=c("cog.id", "year"))

## Pooled Models (Columns 1 and 2)
m1 <- plm(state.rev.pp ~ state.lobby.lag1 + log.pop + state + year, data = pdta, model = "pooling")
## Cluster SEs by city
coef1 <- coeftest(m1, vcov = vcovHC(m1, method = "arellano", type = "HC1", cluster = "group"))

m2 <- plm(state.rev.pp ~ state.lobby.lag1 + log.pop + state.year, data = pdta, model = "pooling")
## Cluster SEs by city
coef2 <- coeftest(m2, vcov = vcovHC(m2, method = "arellano", type = "HC1", cluster = "group"))


## City FE Models (model = "within") (Columns 3 and 4)
m3 <- plm(state.rev.pp ~ state.lobby.lag1 + log.pop + year, data = pdta, model = "within")
## Cluster SEs by city
coef3 <- coeftest(m3, vcov = vcovHC(m3, method = "arellano", type = "HC1", cluster = "group"))

m4 <- plm(state.rev.pp ~ state.lobby.lag1 + log.pop + state.year, data = pdta, model = "within")
## Cluster SEs by city
coef4 <- coeftest(m4, vcov = vcovHC(m4, method = "arellano", 
                                    type = "HC1", cluster = "group"))


## Get DV mean
transfer.mean <- paste0("$", round(mean(pdta$state.rev.pp)))


## Table 1: City Lobbying and State Revenue Transfers Across the U.S., 2007 and 2012
stargazer(coef1, coef2, coef3, coef4, 
          star.cutoffs = c(0.05), type = "text",
          title = "City Lobbying and State Revenue Transfers Across the U.S., 2007 and 2012.", 
          dep.var.caption= ("State Revenue Transfer to Cities (Per Capita), $t$"), 
          keep = c("state.lobby.lag1", "log.pop"),
          covariate.labels = c("Lobby, $t-1$", "Population (Log), $t$"),
          add.lines = list(c("State FE", "Y", "", "", ""),
                           c("Year FE", "Y", "", "Y", ""),
                           c("State x Year FE", "", "Y", "", "Y"),
                           c("City FE", "", "", "Y", "Y"),
                           c("Observations", N(m1), N(m2), N(m3), N(m4)),
                           c("\\# Cities", length(levels(pdta$cog.id)), length(levels(pdta$cog.id)),
                             clusters(m3), clusters(m4)),
                           c("State Transfer Mean", rep(transfer.mean, 4))),
          notes.append = FALSE, notes.align = "l", notes.label  = "",
          notes = "Robust standard errors clustered by city. $^{*}$p$<$0.05")
